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METHOD FOR IMPROVING DEPH DISCRIMINATION IN 



OPTICAL REPRODUCTION SYSTEMS 



[0001] The invention is directed to a method for depth discrimination in optically imaging 
systems. In particular, it is applicable in light microscopy for achieving an improvement in 
image quality in the examination of three-dimensionally extending objects. 

[0002] A problem arising in conventional light microscopy when examining three- 
dimensionally extending objects, i.e., objects whose extension along the optical axis is 
greater than the depth of focus of the objectives that are used, consists in that the in-focus 
image from the region of focus is superposed with extra- focal image components which are 
accordingly imaged out of focus. 

[0003] In order to overcome this problem, it is known to use confocal imaging in which 
the light originating from the area outside the region of focus is excluded by means of a 
pinhole and therefore does not contribute to the image. In this way, an optical section is 
formed. This confocal point-type imaging requires that the object is scanned in the image 
plane in order to obtain an image. This scanning can be carried out either by means of 
scanners in laser scanning microscopes or by means of Nipkow disks. 

[0004] By recording a plurality of optical section images in different focus positions, a z- 
stack can be obtained so that the object may be displayed three-dimensionally. 

[0005] Another method for generating optical sections is to use structured illumination. 
This was first stated by Meir Ben-Levy and Eyal Pelac in WO 97/6509. This method was 
improved and expanded as described by Tony Wilson et al. in WO 98/45745 and by Volker 
Gerstner et al. in WO 02/12945. The disclosure of these three publications is expressly 
referenced herein. v 

[0006] In WO 97/6509, the object is illuminated by a periodic structure (sine grating or 
rectangular grating) and an image of the object is recorded by a camera, digitized and stored 
in a storage. The periodic structure is subsequently displaced within the image plane in such 
a way that the phase position of the structure is changed and an image is again recorded and 
stored. This process (displacement, image recording, storage) can be repeated many times. 



A section image is then generated by calculating from the existing images. The indicated 
mathematical formulation is a Fourier expansion, which leads to a complicated formula 
apparatus. 

[0007] WO 98/45745 shows a simpler formula for the section images which can be 
derived by simplifying the formula from WO 97/6509 in case of equal phase displacements 
between the individual recordings. 

[0008] Realizing the physical boundary conditions required for the application of the 
indicated method proves very difficult in practice. For example, variation in lamp brightness 
between the different recordings leads to stripe artifacts in the generated section images. 
With fluorescing objects, additional problems occur due to the time-dependent fading of the 
fluorescent dyes, which likewise results in errors. The necessary constancy of the individual 
phase displacement steps cannot be maintained in practice. 

[0009] Therefore, it was suggested in WO 02/12945 to compensate for the influence of 
lamp brightness that varies over time by coupling out part of the light serving to illuminate 
the object, registering the intensity and subsequently scaling the individual recordings. In 
order to take into account unequal phase displacement steps, a system of equations (Equation 
22, op. cit.) is indicated. To compensate for fading or bleaching, instead of the necessary 
minimum of three recordings per section image, it is suggested that six recordings are 
registered in the sequence 1-2-3-3-2-1, two section images are calculated (from 1 - 2 - 
3 and 3-2-1) and the average is determined therefrom. 

[0010] A considerable expenditure on instrumentation is required in order to realize these 
suggestions. Further, the recording of additional images prolongs the required recording time 
and therefore also increases loading of the sample by illumination with the fluorescence 
excitation light. 

[001 1] It is the object of the invention to overcome the disadvantages of the prior art and 
to provide an improved method for determining optical section images. This object is met by 
a method for improving the depth discrimination of optically imaging systems according to 
the main claim. Advantageous developments are indicated in the subclaims. 



[0012] The invention will be described more fully in the following with reference to the 
drawings. 

[0013] Fig. 1 shows the basic construction of a microscope with structured illumination; 

[0014] Fig. 2 is a schematic view of the bandpass filter, according to the invention, in the 
Fourier plane; 

[001 5] Fig. 3 is a schematic view of the position of the images of the periodic structure; 

[0016] Fig. 4 shows a flowchart for determining correction values through linear 
optimization. 

[0017] Fig. 1 is a simplified optical schematic illustrating the structured illumination with 
transmitted illumination by way of example. The imaging beam path (image-forming beam 
path) is shown. A one-dimensional periodic structure (transmission grating) (3) which is 
located in a focus plane of the optical arrangement shown in the drawing is illuminated by a 
light source (1) through collector optics (2). The grating is followed in light direction by a 
plane-parallel glass plate (4). The angle of the plane-parallel plate relative to the optical axis 
can be adjusted in a defined manner. The structure is imaged in the specimen plane (7) by 
the illumination-side optics (5, 6) (condenser). Imaging is effected by the light coming from 
the specimen through a pair of lenses (8, 9) (objective and tube lens) in the focus plane (10) 
following the latter in which, e.g., a CCD matrix of a digital camera can be arranged. By 
tilting it in a defined manner, the plane-parallel glass plate (4) serves to displace the image of 
the grating structure (3) on the object in the specimen plane (7). For incident fluorescence 
observation, the objective (8) preferably serves at the same time as a condenser. As was 
described in WO 02/12945, the resulting brightness distribution in the object is registered in 
at least 3 positions of the glass plate (4) by means of the digital camera. In this connection, 
the brightness distribution l { (x, y) for a sine/cosine grating can be stated in a simplified 
manner by the following formula: 

I t (x, y) = I 0 (x, y) • (1 + m(x, y) - cos(^ (x, y) + $ (x, y)) , (1) 

where i=0 . . . N-l is the ith phase position of the projected grating, and N is the quantity of 
recordings, m(x, y) is the modulation depth of the object (and therefore the image 
information sought at point x, y), and (/> \ are the phase values. This equation contains the 



three unknown variables Io, m and (p o. Accordingly, these variables can be determined by 
means of at least three measurements with specifically varied (/> \ (i = 1, 2, 3). The solution 
can be found from the measurements through a least square formulation. For this purpose, it 
is possible to represent Equation (1) in a more compact form and to rewrite the cosine 
function using the addition theorem: 

/(x, y 9 fa ) = a 0 (x, y) + a x (x, y) - J[ (fa ) + a 2 (x, y) • f 2 (fa ) , (2) 
where 

Mfa) = cosfa 
f x (fa) = smfa 
a Q (x,y)^I 0 (x,y) 

a x (x, y) = I 0 (x, y) • m(x, y) • cos <f> 0 (x, y) 

<* 2 (x, y) = -/<> (*> y) * m ( x > ' sin (*> y) 

[0018] The functions fi and f2 therefore depend only on the phase displacements $ , which 
can be freely selected in principle. For a structure that is not sinusoidal but is periodic, the 
brightness distribution Ii (x, y) can also be approximated by series expansion. The principle 
of calculation remains basically the same. Expressed in matrix form, the least square solution 
is: 



M-a=b , 
where 

f 



M - 



N N 
N N N 

K N N N 



(3) 



(4) 



where N is the quantity of measurements (in this case, phase steps) and 

(5) 



a = 



r ao(.x,y)> 
a\{x,y) 
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and 



f 



X/(*,M) 



N 



b = 



£[/(*,M)/i(4)] 



(6) 



N 



£[/(*,M)-/ 2 (A)] 



V N 



J 



[0019] The solution to this matrix equation is obtained mathematically by inverting the 



[0020] Accordingly, the unknown variables being sought can be determined from the 
measurements, represented by vector b, and the phase displacements of the grating. In 
particular, the modulation depth at point (x, y) is given by: 



[0021] Since the grating contrast and, therefore, by necessity, the modulation depth in the 
object decreases with defocusing of the object, Equation (8) is a conditional equation for 
obtaining a depth-discriminated image (optical section). Equations (2) and (3) provide a 
universal approach that allows the phase positions $ to be adapted to the actual conditions. 
Accordingly, the number of steps and the positions in the object space (location of the 
imaged grating lines) can be freely selected. 

[0022] Each optical system, according to its optical transfer function, transfers spatial 
frequencies up to a threshold value that is determined by the numerical aperture NA of the 
optics and the light wavelength X,: 



matrix M. 



<*o(x 9 yy 

a x {x,y) =M~ l 'b 



(7) 




(8) 



_ AnNA 

max ~ -5 



[0023] As a result, the higher harmonics of the grating period of a grating which is 
imaged in the object space and which is not a pure sine grating are transmitted up to this 
threshold frequency k max and are found in the Fourier-transformed image of the object as 
local maxima. These components of the grating are also contained in a section image 
determined according to the plan detailed above and result in stripe artifacts. These artifacts 
can be eliminated in that the harmonics of the projected grating are localized in the spatial 

£ , ,-, .-v r> 4 i •. Kir I7sviivi or- frnncfrvrmoti r»n onrl arp rf*mr\\TF*c\ \\\r a cnf*r»ifif Cf^ttin cr r\€ fl 1 tpr<v 

11 CV^UV/ilW^ 0|J^/WIA U-111 Ojr A V/Vili^/A uuiiuiwiniunwii uxjivi *- * -«-" » * J *- -'^ ~ ~~ ^ _ „„. ^ 

(bandpass filters, bandstops). A subsequent inverse transformation then results in images that 
are free of artifacts. This method is shown schematically in Fig. 2. In the Fourier plane (1 1), 
the higher harmonics of the grating structure (3) are designated by (12), the zero point is 
designated by (13) and corresponds to the DC component of the illumination, i.e., a uniform, 
unstructured illumination. By introducing bandpass filters corresponding to the harmonics 
(12) into the beam path, these areas do not contribute to the imaging and are therefore not 
registered by the digital camera. 

[0024] With respect to fluorescing objects, there is another reason for the occurrence of 
stripe artifacts, namely, fluorescence bleaching (photobleaching, fluorescence fading). In this 
case, the imaging of the structure (e.g., of the grating) in the object results in regions of 
differing degrees of bleaching of the fluorescence and therefore, finally, in stripe artifacts in 
every image recorded subsequently. The frequency components of the artifacts in the 
spectrum of an image of this type are particularly strong in the fundamental frequency of the 
projected structure. By removing the spatial frequency components of the projected structure 
in the Fourier-transformed image by means of corresponding bandpass filters and subsequent 
inverse transformation, an artifact-free image can also be generated in this case in principle. 

[0025] Prevention of these artifacts occurring in fluorescence illumination can also be 
achieved by determining and making use of the bleaching characteristics of fluorescing 
objects in the sequential recording of images of the same object. In this case, the illumination 
results in bleaching of the object at the illuminated points for every image recording. Apart 
from the illumination intensity, the degree of bleaching also depends on the photochemical 
and photophysical characteristics of the fluorescing object. By determining the fluorescence 
intensity of the object through integration of the intensities over all object points or selected 
object points (x, y), it is possible to rescale the image sequence after the images have been 
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recorded. For this purpose, the intensities at the image points (x, y) are scaled to a maximum 
value either by means of an analog camera and analog-to-digital converter or a CCD camera 
using image information stored in a computer by determining the quotients of the 
fluorescence intensities of the images by pairs. In this way, the fluorescence loss due to 
fluorescence bleaching is taken into account. In addition, according to Fig. 3, regions of 
interest (ROI) can be defined which are equally illuminated in both image recordings and 
which subsequently fluoresce; the degree of the reduction in fluorescence in these regions is a 
measurement for the bleaching rate and can therefore be used for scaling the different 
recordings of the brightness distribution. These ROIs can be defined by the user by means of 
corresponding input devices or can be defined automatically (from the respective position of 
the imaged structure). In Fig. 3, the ratios in a first recording of the brightness distribution on 
the object is shown at (14) and the ratios in a second recording with changed phase position 
of the imaged structure (16) are shown at (15). The regions of interest (17) are 
advantageously defined in such a way that uniform illumination is carried out in both 
recordings. 

[0026] The method according to the invention allows optical section images to be 
determined practically for any selected phase displacements $ . In order to prevent errors, the 
phase displacements^, must be determined with high accuracy. This can be ensured in a 
simple manner by the calibration method described in the following, in which the influences 
of the entire system are taken into account. 

[0027] For this purpose, a reference image (object with superimposed grating) is initially 
recorded and registered. A slight phase displacement of the grating to an initially unknown 
phase position of the grating is then carried out (this "displacement per control signal" is to 
be calibrated first) and a new image recording is made. The two images that are obtained are 
then compared by computation. Without limiting generality, this can be carried out by 
differentiating, summing, or in any other way that results in a merit function. The initial 
result is generally a striped image. The steps comprising displacement of the grating by a 
small amount, image recording, and comparison to the reference image are repeated until the 
merit function reaches the extreme value, that is, e.g., a differential image comprises only 
background noise or a sum image reaches maximum values. When the latter is 
accomplished, the calibration is terminated and a control value obtained in this way can be 



kept in a storage medium so that it can be retrieved at a later time. Therefore, this value 
corresponds exactly to the displacement of the grating by a full period. Alternatively, it is 
also possible to assess a sum image for uniform brightness distribution (i.e., disappearance of 
stripe structure). In this case, this value corresponds to the displacement of the grating by 
half a period. This procedure is preferably carried out with a mirror as object. The described 
process can be carried out manually or automatically. 

[0028] In the following, an alternative method for determining the influences of variations 
in lamp brightness, bleaching of the object in fluorescence illumination, and non-sinusoidal 
gratings is described with reference to Fig. 4. 

[0029] For this purpose, the change in the illumination over time will be described by 
intensity-proportional scalar factors 0, > 0 which are taken into account in every successive 
recording (the dependency of the variables gi, o i9 etc. on location (x, y) is no longer explicitly 
indicated in the following for the sake of simplicity): 



[0030] The modeled observation g, is the product of the ideal, consistent phase image o,\ 
In order to describe a faulty system completely, only N-l factors of 0% need to be determined, 

[003 1] The bleaching of the object in fluorescence illumination depends on the effective 
radiation intensity over time. This varies spatially because of the projected grating and 
therefore must be handled in such a way that the grating function is taken into account. A 
simple description is possible when the bleaching is assumed to be a linear function of the 
illumination intensity: 



[0032] The vector 0 < < 1 comprising spatially varying factors describes the attenuation 
as a function of the shape of the irradiation intensity and a selectable measurement 
1 > d > 0 for the bleaching when multiplied component-wise by the observations o,- that have 
not yet undergone bleaching: 



gi = Opt 



(9) 



k f =(l-rfo M k H )^ 0 k y ,/>^ , where k 0 =1 



(10) 



& = o/kf 



(11) 



[0033] It follows from (10) that every observation is always multiplied by the factors of 
all preceding phase recordings; the latter are already bleached in part and overlap with the 
instantaneous phase. Alternatively, k/ can also be described exponentially, i.e., 
corresponding to the decay curve of a determined fluorochrome. 

[0034] Non-sinusoidal grating shapes cause stripe artifacts particularly for integral 
multiples of the grating frequency. With exact knowledge of frequency, phase position and 
curve shape of the grating in the individual recordings, it is possible to generate a sinusoidal 
shape. A correction vector \corr similar to (10), but in this case along the grating periodicity, 
is considered in (12). 

Icorr = , /Si " - (12) 
J grating s 

[0035] The simulation of the observed grating function is characterized by f grating- An 
ideally sinusoidal intensity curve is expressed by f sin . The constant s in the denominator 
should take over the division by very small values of f grating- The simulation of the observed 
grating function can be carried out in different ways. A simple possibility for approximation 
is provided, e.g., through the use of a Fourier series for the trapezoid function. 

fgratingM) = - + — ^ ^ ' sin(2*l-x) (13) 

[0036] The edge steepness of the trapezoid transitions can be varied by 6, x is the spatial 
coordinate, and M is the quantity of Fourier coefficients used. Alternatively, more exact 
models can also be used. They require convolution of the true grating function with a lateral 
component of the PSF of the microscope. A phase-synchronous synthesizing of the 
correction vector Icorr is crucial. When this is ensured, a corrected result for the observation 
can be achieved similar to (11) by a simple component- wise multiplication. 

giCorr = gUiCorr (14) 

[0037] The corrected phase images are now obtained from: 

J. = 0/ = Zfocorr ^ where 0. > Qnk . > 0 (15) 
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By inserting (15) in Equations (3) or (8), a general formula is obtained for reconstruction 
which, moreover, contains the parameters ft, d and b. 

[0038] The stated problem can be solved, for example, by converting to an extremal 
problem. The parameters ft, d and b can be determined with the aid of numeric optimization 
using a merit function. 

M{9 t \d\b) = oc 0 \F{a} 0 f +a, \f {a)J\ ^a 2 \F{a}J^^^aJ\F{a}J^min(i6) 

In this case, F{a} Q} is a component of the Fourier transformation of the resultant vectors . 

from Equation (5). Functional transformations other than the Fourier transformation are also 
suitable in this connection. 

[0039] The merit function (Equation (16)) is mentioned here only by way of example. 
Other merit functions are also applicable. The aim is to vary ft, d and b in such a way that 
(16) is minimized. Many methods are available for this purpose which are by no means 
limited to gradient methods or line-search methods. Further, every transformation coefficient 
is weighted by a coefficient a,. Accordingly, the algorithm can be adapted to different signal- 
to-noise distances or preferred frequencies. An advantageous value in this respect for 
weighting the DC component is given by the following equation: 

\ F{5}i \ 2 + \F{a} 2 \ 2 + ... + \F{a} n \ 2 
a 0 oc ■ ! (17) 

l F ( 5 } 0 | 

[0040] The general operation of the correction apparatus described herein is illustrated by 
the schematic flowchart in Fig. 4. The use of projections is advantageous for efficient 
determination of parameters ft, d and b. However, the determination of parameters is not 
limited to this one-dimensional case, but functions in the plane in an analogous manner. 

[0041] With a fixed focus position of the camera, at least three brightness distributions 
(18), (19), (20), and possibly more (21), are registered in different phase positions i=l, 2, 3 
. . ., N of the grating structure imaged on the object. These two-dimensional brightness 
distributions are converted, respectively, to a one-dimensional distribution (22), (23), (24), 
and possibly (25), by means of projection. With initial values for the parameters ft, d and b, 



corrected phase images are calculated according to formula (15) and a first image of the 
optical section is calculated (26) therefrom by means of the formulation from formula (3). 
The initial values of parameters 0;, d and b can be entered by the user or also estimated 
automatically. Known minimizing or maximizing tools of a merit function can be used for 
this estimation (e.g., Simplex, Marquardt, or the like). The merit function according to 
formula (16) is determined from this result in (27) and a test is made in (28) to determine 
whether or not a minimum has been reached. If this minimum has not yet been reached, a 
new estimation of the parameters d and b is carried out (29) and continued with (26). 
When the minimum of the merit function has been reached, i.e., the value of parameters 9 h d 
and b can be increased again by further variation of parameters 0„ d and b, the optimal 
parameters are determined and the optical section image is determined in (30) by means of 
formulas (15), (3) and (8) and outputted in (31). Subsequently, the focus position can be 
changed and another optical section image can be determined, and so on, by means of the 
described method. By means of this scanning of the object in z-direction, a z-stack of section 
images and, accordingly, a three-dimensional representation of the object can be generated. 

[0042] Another approach for taking into account bleaching effects in the object in 
fluorescence illumination is based on determination of a local correction function. For this 
purpose, a correction function Ki(x) is assumed: 

whose local value is formed by a grating period r = IkIco. The formulation of this function 
can also take into account the dependency on y, although this is omitted in the present case 
for the sake of simplicity. Using a spatially variable bleaching function Ofa) and by suitable 
conversions with the approximation (see also formulas (1) and (9) 
g.(;c) = ^.(x)/ 0 (x)[l + m cos{o)x + fy)] , it can be shown by the following equation 

f + r/ ; ] i 0 (t)dz + m r r/ 2 1 0 co« + n )d$ 

fC.(x) — j x-t/2 Jx-t/2 

that this correction function is proportional to \/9,{x) with uniform bleaching function. 
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[0043] It has proven advantageous to average over more than one grating period. Further, 
it has been shown that this correction function tends to spike in the neighborhood of edges in 
the bleaching function and can therefore generate additional errors in the results. In order to 
remedy this, it has proven successful to detect the occurrence of spiking by comparing to a 
mean value and to use the global value 1/ 0t previously determined from formulas (9) and 
(16) at these points instead of the falsified value. The corresponding formula is then given by 



[0044] For this purpose, the threshold 8 is advantageously set as a percentage of the 
variation of the local variability of the correction function within the image, e.g., 5%. Instead 
of 1/0,-, interpolated values of k;(x) which were calculated below the threshold in the two- 
dimensional space can also be used. 

[0045] This method accordingly takes into account bleaching phenomena which vary in a 
spatially-dependent manner, e.g., due to different fluorophores or varying characteristics 
thereof 

[0046] The realization of the invention is not limited to the embodiment examples shown 
herein. Further developments by persons skilled in the art or other approaches for taking into 
account the brightness variations of the light source or of the bleaching of the object in 
fluorescence illumination remain within the field of the invention. 
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